import sys

# argv: <ref_fasta> <alt_fastq> <svm_weight> <output> <threshold>

def rev_com(seq):
    rev = {}
    rev['A'] = 'T'
    rev['T'] = 'A'
    rev['C'] = 'G'
    rev['G'] = 'C'
    rev['N'] = 'N'
    seq_rev = ''
    for i in range(len(seq)):
        seq_rev += rev[seq[len(seq) - 1 - i]] 
    return seq_rev

kmer_len = 10
tf_n = 94

out = open(sys.argv[4],'w')

# Reading thresholds
tf_name = [] 
tf_thresh = []
th = open(sys.argv[5], 'r')
line = th.readline()
while line:
    tf_name.append(line.strip().split('\t')[0])
    tf_thresh.append(float(line.strip().split('\t')[1]))
    line = th.readline()
th.close()

# scanning possible k-mer
ref = open(sys.argv[1], 'r')
alt = open(sys.argv[2], 'r')
line_ref = ref.readline()
line_alt = alt.readline()
kmers = set()
while line_ref and line_alt:
    seq_ref = line_ref.strip().upper()
    seq_alt = line_alt.strip().upper()
    if seq_ref[0] == '>':
        line_ref = ref.readline()
        line_alt = alt.readline()
        continue
    for i in range(31):
        kmers.add(seq_ref[i:i+10])
        kmers.add(seq_alt[i:i+10])
    line_ref = ref.readline()
    line_alt = alt.readline()
ref.close()
alt.close()
print('kmer num : ', len(kmers))

# reading svm weights
svm_m = open(sys.argv[3], 'r')
svm_score = {}
line = svm_m.readline()
while line:
    kmer = line.strip().split('\t')[0]
    score = line.strip().split('\t')[1:]
    kmer_rev = rev_com(kmer)
    if kmer in kmers:
        svm_score[kmer] = score
    if kmer_rev in kmers:
        svm_score[kmer_rev] = score
    line = svm_m.readline()
svm_m.close()

# scanning sequence
ref = open(sys.argv[1], 'r')
alt = open(sys.argv[2], 'r')
line_ref = ref.readline()
line_alt = alt.readline()
while line_ref and line_alt:
    seqid = line_ref.strip()[1:]
    line_ref = ref.readline()
    line_alt = alt.readline()
    seq_ref = line_ref.strip()
    seq_alt = line_alt.strip()
    deltasvm_score = [0 for i in range(len(tf_name))]
    for i in range(tf_n):
        for j in range(31):
            if seq_ref[j:j+10] in svm_score:
                deltasvm_score[i] -= float(svm_score[seq_ref[j:j+10]][i])
        for j in range(31):
            if seq_alt[j:j+10] in svm_score:
                deltasvm_score[i] += float(svm_score[seq_alt[j:j+10]][i])
        pred = "None"
        if abs(deltasvm_score[i]) > tf_thresh[i]:
            pred = "Gain" if deltasvm_score[i] > 0 else "Loss"
        out.write('{0}\t{1:.6f}\t{2}\t{3}\n'.format(seqid, deltasvm_score[i], tf_name[i], pred))
    line_ref = ref.readline()
    line_alt = alt.readline()
ref.close()
alt.close()
